Frequency-dependent spin susceptibility in the two-dimensional Hubbard model 
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A Quantum Monte Carlo calculation of dynamical spin susceptibility in the half-filled 2D Hub- 
bard model is presented for temperature T = 0.2t and an intermediate on-site repulsion U — At. 
Using the singular value decomposition technique we succeed in analytically continuing the Mat- 
subara Green's function into the real frequency domain and in deriving the spectral representation 
for the longitudinal and transverse spin susceptibility. The simulation results, while contradicting 
the random-phase approximation prediction of antiferromagnetic long-range order at this temper- 
ature, are in agreement with an extension of a self-consistent renormalization approach of Moriya. 
The static susceptibility calculated using this technique is qualitatively consistent with the lj — > 
simulation results. 

PACS numbers: 74.20.-z, 74.20. Mn, 74.25. Dw 



The dynamical magnetic susceptibility in strongly cor- 
related electron systems (SCES) has been at the cen- 
tre of attempts to explain the pairing mechanism in 
high temperature superconductivity (HTSC) |^ as 
well as to understand the normal-state properties of the 
HTSC compounds The frequency dependence of 

the magnetic susceptibility can be measured in neutron 
scattering experiments, thus allowing the assessment of 
the relevance of various low-energy models of SCES to 
HTSC, and the validity of their approximate solutions 
in various temperature and doping regimes. Quantum 
Monte Carlo (QMC) simulations, among other numeri- 
cal techniques, provide an alternative test of these models 
and of the applicability of the analytical approximations. 
These simulations supply direct information only about 
the imaginary-time dependence of the correlation func- 
tions. Therefore, most of the information derived from 
them has been concerned with the static susceptibility 
1^. Analytic continuation of QMC data into the real- 
time domain has been performed recently to obtain the 
single-particle spectral weight function and two- 
particle 1^ Green's functions. However, a comprehen- 
sive dynamical description of the SCES, particularly of 
their collective excitations, based on QMC simulations, 
remains a largely unexplored area. In this letter we re- 
port results of analytic continuation by means of singu- 
lar value decomposition (SVD) to evaluate the dynami- 
cal spin susceptibility. These results allow us to examine 
spectral characteristics predicted in the spin-density wave 
(SDW) treatment for a 2D Hubbard model at finite 
temperature. 

Interpretations of QMC results for magnetic suscepti- 
bility have been mostly related to random-phase approxi- 
mation (RPA) calculations which lead to an SDW ground 
state Comparisons have been made between finite- 
lattice simulations and RPA results for an infinite system 
thus neglecting the finite-size effects. These latter 
calculations overestimate the value of the Neel temper- 



ature which would be only enhanced in a calcula- 
tion done on a finite lattice Even without this en- 
hancement, the fit of the QMC data to the RPA predic- 
tion requires an artificial correction to the bare Coulomb 
repulsion in the RPA expression for the susceptibility 
§ , ||ll[] , (l2| . We explored, therefore, an application of 
an extended version of the self-consistent renormalization 
(SCR) theory which is known to provide a more accurate 
description of weak antiferromagnetism [p^ , and found 
that for a finite system this theory agrees qualitatively 
with the QMC results with no parameter-adjustment. 

We apply a finite-temperature QMC algorithm simu- 
lating the partition function of the Hubbard model as a 
path integral of a euclidcan field theory. The temperature 
is represented by an additional compact dimension whose 
extent is equal to the inverse temperature (3. Matsubara 
Green's functions Gab{t) for operators A and B, satis- 
fying for imaginary-time separation t G(r + 0) = ±G{t) 
are evaluated as ensemble averages of the type: 



Gab(t) 



J[D^]A{t)B{0) exp(5[$]) 
cxp(S'[$]) ' 



(0 < r < /3) 



(1) 



where $ are the generalised coordinates of the action S. 
Analytic continuation of these functions into real time 
amounts to a solution of the following ill-posed inverse 
problem: 



G(r) 



duj 



1 ± e-'^f^ 



X"iuj); (0<r</3). (2) 



In what follows x"{^) is the imaginary part of the com- 
mutator retarded Green's function, which determines the 
sign in the denominator (— ) and the periodicity of the 
boundary condition. The solution for x"i'-^) '^^s found 
using the singular value decomposition technique de- 
scribed in Ref. |^. 
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We calculate spin-correlation functions derived from 
QMC simulation data for Matsubara Green's functions 
of the type: 

Gs^S^ (fc, r) = {{n^ (fc, r) -n^ (fc, r))(nt (-fc, 0)-ni{-k, 0))) 
Gs+s- (k, t) = (c| (fc, T)q (fc, r) cj (-fc, 0)ct (-fc, 0)) (3) 

where cj. (c^) is an electron creation (annihilation) op- 
erator with spin projection cr, Ua = c\ca and k is the 
lattice momentum. The simulations were performed for 
a single-band 2D Hubbard model with on-site Coulomb 
repulsion U = At {t being the nearest-neighbour hopping 
parameter) at simulation temperature 1/(3 = 0.2t and 
lattice-sizes 8^ and 10^. The models were simulated at 
half- filling. We solve Eq.(|2|) with its left-hand-side given 
by Gs's^{k,T) and 6*5+5- (fc,T), its solution x"{k,i^) 
being the imaginary part of the dynamical longitudinal 
(x^^(fc, Lu j) and transverse (x^ {k, uj)) susceptibilities re- 
spectively. 

In Fig.l we present the results for x"(fc, w) for various 
values of the lattice momentum k. The data is based 
on 4200 measurements on a 8^ lattice and 4100 measure- 
ments on a 10^ lattice. The ensemble averages and the 
statistical errors were estimated using bootstrap analy- 
sis. There are no qualitative differences due to finite-size 
effects between the two sets of results, which suggests 
the possibility of extrapolating the essential behaviour 
to the thermodynamical limit at this temperature. The 
analytic continuation has been based on reconstruction 
with seven singular functions of the kernel in Eq. (^) . This 
number corresponds to seven singular values whose ratio 
to the largest one is not smaller than the average er- 
ror (0(10"^)). The Matsubara Green's functions were 
evaluated accordingly at seven values of r, spaced corre- 
spondinly to the zeros of the seventh singular vector. 

A sum-rule used in Ref. for the t — J model can 
be applied in our case to examine the extent of double- 
occupancy in the given temperature and coupling regime. 
With our normalization (see Eq.(0)) this sum-rule reads: 



l^°°d^X^coth^x'"^(fc,c.) 

k 



{St 



(4) 



where N is the number of the lattice sites. Our results for 
the average on-site spin projection give ((,Sf)^^ w 0.74 

which suggests that 13% of the sites are doubly occupied. 
This violates significantly the constraint assumed in the 
t-J model and undermines the validity of extrapolating 
results derived for the t-J model to the Hubbard model 
in this regime of coupling and temperature |lj] . 

To compare our results with those derived from the 
SDW-RPA treatment we examine our results with re- 
spect to the evidence for long-range antiferromagnetic 
ordering. The SDW-RPA predicts, for the size of the sys- 
tems simulated, a Neel temperature T^^^ ~ 0.7t, which 



is above our simulation temperature. The same treat- 
ment predicts also an SDW gap A = 1.38t at our simula- 
tion temperature. This gap would result in the vanishing 
of the time-ordered Green's function at the nesting vector 
g=(7r,7r): x""" ^(Q, ^) for a; < 2A §. 

We find that x"(fc, t^) for small values of uj has a clear 
peak at fc = Q which leads also to a sharply enhanced 
real part of static susceptibility x'iQi'^ = 0) calculated 
using the Kramers-Kronig relations (Fig. 2). However our 
results contradict two essential SDW predictions: the 
spontaneous breaking of rotational symmetry and the ex- 
istence of a gap in the spectral function for the magnetic 
correlations. The approximate preservation of the rela- 
tion x^^(fc, w) = 2x"'"~(fc, w) (Fig.l) indicates that the 
rotational symmetry is not violated, while the results for 
^i/zz T^Q^ j^-j (pjg 3) show no evidence for the SDW gap. 
This latter discrepancy with the SDW prediction cannot 
be explained by finite-size effects, since the sign of the 
gap is invariant under the change of the staggered mag- 
netization. Therefore if antiferromagnetic order existed 
in the thermodynamical limit we would expect a mini- 
mum (if not vanishing) of x"^^'^{Q7^) for small values 
of Lu. This is clearly absent in our results. 

We observe therefore that our results cannot be ad- 
equately explained within the SDW-RPA scheme. The 
standard RPA result yields for T > T^ccV 



Xi^pa(A:,w) 



1 



-C/xJ-(fc,c^)' 



(5) 



We note that the non-interacting susceptibility Xo(QiO) 
calculated for a two-dimensional finite system at half- 
filling exhibits a 1/T divergence as T — > 0. Therefore, 
there will always be a pole of the susceptibility at some 
nonzero temperature T^^^, indicated by Uxo{Q, 0) = 1, 
with a subsequent divergence of the total free energy of 
a finite system. This unphysical result is one of the 
consequences of an inherent shortcoming of the SDW- 
RPA treatment, namely, its failure to represent in a self- 
consistent manner the contribution of the spin fluctua- 
tions to the free energy. The free energy can be expressed 
as a functional in the dynamical susceptibility [p^ (see 
Eqs.(|),(|)). On substituting expression (||) into this 
functional and subsequently differentiating it twice with 
respect to the staggered magnetization one fails to re- 
cover the static limit of (|^). This inconsistency has been 
pointed out |l^ as the reason for predicting too high val- 
ues for the Neel temperature of weak-antiferromagnetic 
metals. These arguments suggest the need to explore an 
analytic approach which would more adequately describe 
the thermodynamics of the system, if we are to try to re- 
late our QMC results to any analytical concepts. 

Such an approach has been developed as the theory of 
SCR by Moriya et al. ||l^ . A renormalization function A 
is introduced which relates the exact sum of two-particle 
irreducible diagrams Xum i^, 1^) to the free susceptibility 
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Xom(^''^) ^ given background staggered magnetiza- 
tion M: 



1 + \uM{k,uj) 



(6) 



This renormalization function enters the calculation of 
the full susceptibility: 



X^Mik. c^) A(fc + Q) + Ujxtu jk, k + Q, iu)]^ 
A{k)A{k + Q) - [Ux^j^fik, k + Q, 

A{k) = l + XuMik,iu)-Ux^Mik,^)- (7) 

Here the free umklapp susceptibility at staggered mag- 
netization M is denoted by Xo^(fc, ^ + Q, w). The free 
energy of the system with N^i electrons is given by: 

F{U, M) = Fo{M) + - - M^j + aF{U, M) (8) 

where the first two terms represent mean-field free energy 
and the third one is the correction due to the collective 
excitations which can be expressed via Xu~\i{k,ijj) as: 



aF{U,M)-^ 



~T dU' 
Jo 



k 



(9) 

where = 27rTz^. The sum in (|^) is dominated by 
Xum{QtO)j which justifies the long-wave static approxi- 
mation Xi/M{k,Lu) — Xum{Q,0)- In this approximation, 
Eqs.(0)-(^ together with the thermodynamical relation 
between the static susceptibility and the free energy: 



XlrMiQ,0) 



1 d^FjU, M) 
2 



A/=0 



(10) 



constitute a set of equations for the susceptibility, auto- 
matically maintaining the self-consistency requirement. 
Eq.(|l0|) can be solved by expanding \um{Q,^) in pow- 
ers of M. In the next to leading order expansion: 
Ac/m(Q,0) « Ao(C/) + {1/2)\2{U)M^ we obtaim 

d^AF 



Mu)^xtM{Q,^)^,,,^, 



dAp \m=o' 



(11) 



Since even within this approximation the thermody- 
namics of the spin fluctuations is accounted for self- 
consistently, the qualitative description of the suscepti- 
bility is likely to be more accurate than the one provided 
by SDW-RPA. Apparently due to the finite-size 1/T di- 
vergence of Xoi Moriya's approximation of A independent 
of U and M, used in the infinite-system calculation 
is too crude for this calculation. By substituting the so- 
lution of Eq.(pTl) for A;7m(Q,0) into (||) at w = we 



obtain the value of (fc, 0). The static susceptibil- 
ity (at T > Tncci) is obtained by replacing Xo {k,0) in 
Eq.(|) by Xf^o 0) (Fig.2). The calculation of the direct 
and umklapp free susceptibility entering Eq. (|^) was done 
for lattice sizes 8^ and 10^. Thus the finite-size effects, 
inevitably present in the QMC simulations, have been 
accounted for in the analytical calculation. We note the 
qualitative agreement between the calculated values and 
the simulation results at the nesting vector and the rest 
of the lattice momenta. The agreement is particularly re- 
markable since there is no parameter adjustment in this 
analytical calculation. Due to the discrete spectrum of 
single-particle excitations on a finite lattice, Xoik, lu) is a 
sequence of (5— functions and becomes smooth only for an 
infinite system. An imaginary regularizer (corresponding 
to a finite life-time for a quasiparticle) has thus to be in- 
troduced to allow the calculation of dynamical suscepti- 
bility using Eq.(0). We found the results to be sensitive 
to the value of this regularizer, therefore, a meaningful 
comparison with the QMC data can be done only for 
static susceptibility. 

Making a standard approximation x'(^+ Q,0) = 
x'(Q, 0)/(l -I- £,^(1^) for small deviations q from the nest- 
ing vector, we obtain the correlation length £ — 3.2a and 
^ = 3.6a (a being the lattice spacing) for the 8^ and 10^ 
lattices respectively. The SCR estimates are £ = 3.9a and 
^ — AAa. We emphasize, however, that since the correla- 
tion length calculated is of the order of the lattice linear 
dimension this estimate can vary as lattice size increases. 

To summarise, the SVD technique is applicable to the 
derivation of the dynamical properties of collective exci- 
tations in SCES. The dynamical susceptibility obtained 
from QMC simulations contradicts the qualitative pre- 
dictions of SDW-RPA, in particular there is no evidence 
for SDW gap formation in the longitudinal time-ordered 
correlation function. Therefore it is unlikely that the 
gaps observed for this lattice size and temperature in 
previous QMC simulations of the single-particle density 
of states can be explained within SDW-RPA as was 
indeed pointed out in Rcf.Q. Our results support the 
evidence for the insufficiency of the SDW-RPA picture 
to describe the spectrum of excitations of the Hubbard 
model ijl^ , jl^ . Their qualitative agreement in the static 
limit with the SCR calculation points to the importance 
of including the paramagnon interaction in the descrip- 
tion of the magnetic properties of the Hubbard model. 
This approach has previously been used to obtain on a 
more phenomenological basis, the dynamics of collective 
excitations in HTSC compounds in the normal phase as 
well as to calculate the pairing potential Our work 
outlines a way to examine these properties directly from 
the lattice-field model by means of numerical simulations 
as well as by a power series expansion in the magnetisa- 
tion within the SCR theory. 
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Figure Captions: 

Fig.l x"""{k,oj) (solid hne) and x"^-{k,u)) (dotted hne) 
vs. to for selected values of the lattice momentum given 
in the brackets in each figure; (a) - lattice 8^, (b) - lattice 
102. 



Fig. 2 Static spin susceptibility vs. lattice momentum; 
(a) - lattice 8^, (b) - lattice 10^. The results of the SCR 
calculation are shown by (•) . 

Fig. 3 Timc-ordcrcd Green's function x"^^^{Q,i^)- Note 
the absence of SDW gap. 
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Fig. 2 
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